#!/usr/bin/perl -w
use strict;
use Getopt::Long;

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

my $usage = <<'USAGE';

USAGE:

    filter_fasta_by_regex_list.pl
                --regex "(ENS[A-Z0-9]+)\t"
                --list list_of_matching_codes
				--reverse
				[--histogram matching_histogram_output_file]
				[--verbose]
				[--help]
	
	Iterates through a file and filters out any non-matching lines
	Codes are extracted from each line using the supplied regular
	    expression and compared to the supplied list.
USAGE


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory
my $regex;
my $list;

# optional parameters
my $verbose;
my $histogram;
my $reverse;
my $help;
GetOptions(
			'reverse'			=> \$reverse,
			'list=s'			=> \$list,
			'regex=s'			=> \$regex,
			'histogram=s'		=> \$histogram,
			'verbose'			=> \$verbose,
			'help'				=> \$help);

die $usage if ($help);
die $usage unless ($regex && $list);

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	get list of matching codes

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

open FILE, $list or die "Error!\nCould not open the file [$list]\n$!\n";
my @file_codes = <FILE>;
chomp @file_codes;
close FILE;

# uniquify array of codes and save duplicates
my %codes;
$codes{$_}++ foreach (@file_codes);

my @duplicates;
my $cnt_duplicate = 0;
foreach (keys %codes)
{
	push(@duplicates,  $_) if $codes{$_} > 1;
	$cnt_duplicate += ($codes{$_}-1);
}

# set hash values to zero: this will hold the number of lines matching that ID
$codes{$_} = 0 foreach (keys %codes);


die "Error!\nNo codes found in [$list]\n" unless @file_codes;

my $cnt_line=0;
my $cnt_filtered_lines=0;
while (<>)
{
	$cnt_line++;
	my @codes = /$regex/g;
	
	# Are there any regex matches on this line?
	if ($verbose && !@codes)
	{
		print STDERR	"Warning!\nNo codes matching [$regex] were found on the ".
						"accession line $.\n";

	}
	
	# Are any of the regex matches in the supplied list?
	my $any_matches;
	foreach my $id (@codes)
	{
		if (exists($codes{$id}))
		{
			++$codes{$id};
			$any_matches = 1;
		}
	}
	$any_matches = !$any_matches if ($reverse);
	print, $cnt_filtered_lines++ if $any_matches;
}

print STDERR "\t$cnt_line\tlines in the original file(s)\n";
print STDERR "\t", scalar @file_codes, "\tcodes in the supplied list\n";
print STDERR "\t", scalar keys %codes, "\tunique codes in the list\n"
				if (scalar keys %codes != scalar @file_codes);
print STDERR "\t", scalar @duplicates, "\tcodes with duplicates in the list\n" if @duplicates;
print STDERR "\t\ti.e. ", $cnt_duplicate, "\tduplicates in the list\n" if @duplicates;
print STDERR "\t", scalar keys %codes, "\tcodes matched lines\n";
print STDERR "\t", '-'x6, "\n";
print STDERR "\t$cnt_filtered_lines\tlines matched a supplied code\n";

# print counts of sequences which match > 1 ID
if ($histogram)
{
	open HISTFILE, ">$histogram" or die "Error!\nCould not open the file [$histogram]\n$!\n";
	my @sorted_keys = sort {$codes{$a} <=> $codes{$b} } keys %codes;
	if (@sorted_keys )
	{
		my $padding_spaces = length($sorted_keys[0]) - 2;
		print HISTFILE "CODE", ' ' x $padding_spaces, "\t", "Count of matching lines\n";
		print HISTFILE "----", ' ' x $padding_spaces, "\t", "-----------------------\n";

		foreach (@sorted_keys)
		{
			print HISTFILE "$_\t$codes{$_}\n"
		}
	}
	
	
	if (@duplicates)
	{
		print HISTFILE "\n\n", scalar @duplicates, " duplicate codes in the supplied list:\n";
		print HISTFILE join ("\n",  @duplicates,  "");
	}


}
